Limit cycles, complex Floquet multipliers and intrinsic noise 
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We study the effects of intrinsic noise on chemical reaction systems, which in the deterministic 
limit approach a limit cycle in an oscillatory manner. Previous studies of systems with an oscilla- 
tory approach to a fixed point have shown that the noise can transform the oscillatory decay into 
sustained coherent oscillations with a large amplitude. We show that a similar effect occurs when 
the stable attractors are limit cycles. We compute the correlation functions and spectral properties 
of the fluctuations in suitably co-moving Frenet frames for several model systems including driven 
and coupled Brusselators, and the Willamowski-Rossler system. Analytical results are confirmed 
convincingly in numerical simulations. The effect is quite general, and occurs whenever the Flo- 
quet multipliers governing the stability of the limit cycle are complex, with the amplitude of the 
oscillations increasing as the instability boundary is approached. 

PACS numbers: 05.40.-a,02.50.Ey,05.45.-a 
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I. INTRODUCTION 



The subject of non-linear dynamics, with its wide 
range of tools and techniques, and its classification of 
the diverse types of behavior encountered, has in the last 
20 or 30 years transformed our understandingof many 
models in the physical and biological sciences All 
these systems are subject to random perturbations, but 
the study of the effects that the noise has on a particu- 
lar system, while still very significant [1, 0) Hi) has not 
been nearly so extensive. Frequently the noise is added to 
the deterministic equations in a fairly ad hoc manner to 
obtain stochastic differential equations of the Langevin 
type. What is less often done is to start from a well- 
defined "microscopic" model defined by either a Markov 
chain or a master equation, and to treat the deterministic 
(macroscopic) limit of the model in a unified framework 
which also incorporates the stochastic elements of the 
problem. In this paper we will develop such a treatment 
for a particular class of problems. Conventional tools 
used in deterministic nonlinear dynamics (for example, 
Frenet frames and Floquet analysis) will turn out to also 
have a role to play in the stochastic version of the model. 

The particular class of problems we shall investigate 
will be those which have a deterministic limit which, at 
large times, approaches a limit cycle in an oscillatory 
manner. That is, trajectories spiral into the limit cycle 
at large times. The motivation for studying such sys- 
tems is the widespread interest that there has been in 
the analogous phenomena in systems which approach a 
fixed point in an oscillatory fashion. In this case, the 
effect of noise is, in many cases, to transform the oscil- 
latory decay into a sustained oscillation about the fixed 
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point. In this way the long-time behavior of the sys- 
tem is no longer a fixed point, but consists of stochastic 
oscillations which have a frequency which may be differ- 
ent to that which appears in the oscillatory decay in the 
deterministic version. The possibility of such an effect 
occurring has been discussed for some time [1, 0], but 
it is only in the last few years that a full quantitative 
analysis has been given. The method has been applied 
to the study of stochastic oscillations in predator- prey 
systems [1, Q , epidemiology [Tol . fill IT^ , chemical reac- 
tions in the cell [l^l , auto-catalytic reactions , among 
others. One of the important aspects of these oscillations 
is that they have an amplitude of order c/^fN relative to 
the deterministic trajectory, where N is the size of the 
system (the maximum number of individuals, molecules, 
etc, that may be put into the system) and c is a constant 
which due to a resonance effect is usually quite large. 
This means that even for relatively large values of TV, 
where the oscillations would be expected to be small and 
stochastic effects negligible, the relative amplitude can be 
of order one, and so the fluctuations may dominate the 
dynamics. This effect is usually referred to as stochastic 
amplification, to avoid confusion with the very different 
effect of stochastic resonance 15]. 

The question that will interest us here is: does a sim- 
ilar phenomenon happen in other contexts, in particular 
when the stable state of a deterministic dynamical sys- 
tem is a limit cycle? Much less work has been done for 
this case as compared with the case of a fixed point, yet 
intuitively we would expect a similar effect to occur. In 
fact, the o nly prev ious studies we are aware of are by 
Wiesenfeld [l^, , who investigated the effects of noise 
on the stability of periodic attractors of various dynam- 
ical systems, such as the driven pendulum. He obtained 
analytical and numerical results on the power spectra of 
fluctuations about the limit cycles of such systems, but 
he adopted the approach that we mentioned above: by 
adding noise to the deterministic equations of motion. 
This is acceptable if the noise is external, as he was en- 
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visaging, but if one wishes to understand the possible am- 
phfication of the underlying fluctuations due to intrinsic 
demographic stochasticity, then one needs to begin with 
the discrete dynamics, as we have already emphasized. 

In a recent paper [Tsl , we have investigated a stochas- 
tic model of the well-known Brusselator system, which 
has a limit cycle in the deterministic limit. However, in 
this model the approach to the limit cycle is not oscil- 
latory. As in the case of fixed points and in the appli- 
cations listed above, a precondition for finding sustained 
coherent oscillations is for the stable limit cycles to be 
approached in an oscillatory manner. For a fixed point 
this condition is that the eigenvalues of the stability ma- 
trix about the fixed point are complex. For a limit cycle, 
the analogous condition is that the Floquet multipliers 
for the equations describing the small deviations away 
from the periodic path are complex. Floquet multipliers 
are found to be real in the Brusselator system, as the 
number of degrees of freedom is not large enough to pro- 
duce complex multipliers, and no coherent amplification 
phenomenon is observed. Part of the motivation for the 
work described in was to put the necessary tools in 
place, and to set the stage, for the investigation of model 
systems in which persistent oscillatory behavior about a 
limit cycle is to be expected. 

We begin with a two-dimensional system. If the system 
is autonomous one of the Floquet multipliers will have a 
value of unity, which, as we will see, implies that the re- 
maining Floquet multiplier has to be real. This means 
that complex Floquet multipliers can only be found in 
two-dimensional systems if they are non-autonomous. It 
is natural to achieve this by imposing an external peri- 
odic driving, so as to induce a limit cycle as the steady 
state. In order to make contact with our previous paper 
[m we here first study the Brusselator forced by an ex- 
ternal periodic driving. As it turns, out this system does 
indeed have complex Floquet multipliers for a range of 
possible values for the parameters of the model. We then 
discuss an autonomous system in three dimensions: the 
Willamowski-Rossler model, first introduced to describe 
chemical chaos. Finally, we consider a coupled set of 
two Brusselator systems as a four-dimensional illustra- 
tion. Although we focus on these particular examples 
in the present paper, the formalism we will develop will 
hold in arbitrary dimensions and it will apply whether 
the system is autonomous or non-autonomous. 

The outline of the paper is as follows. We begin in 
Sec. In] with the forced Brusselator. By avoiding the tech- 
nical complexities of working in general dimensions, ap- 
pealing to some of the results used in our previous paper 
on the unforced Brusselator 18], and not having to use 
the Frenet frame in the analysis, we hope to provide a 
gentle introduction to the basic ideas. In Sec. IIIII we 
extend the analysis to the Willamowski-Rossler model 
which introduces some extra features over and above 
those used in Sec. |TT1 and in Sec. llVl we carry out the full 
analysis for a system in an arbitrary number of dimen- 
sions and illustrate its use on the coupled Brusselator. 



We conclude in Sec. |Vl There are three mathematical 
appendices which cover the details of the formalism and 
some aspects of the calculations for the specific models 
considered in the earlier sections. 



II. FORCED BRUSSELATOR 

In this section we will study the Brusselator system, 
subject to an external periodic for cing . An analysis of 
the unforced model can be found in [1^ , and much of the 
formalism remains unchanged. As it turns out, the intro- 
duction of the forcing actually simplifies some aspects of 
the dynamics as discussed below. While we re-iterate the 
main elements of the formalism and of the notation in the 
present paper, our previous work JJj] may be consulted 
for specific details. 



A. Model definitions 

The Brusselator model is a relatively simple chemical 
system, composed of five different reactants [A, B, C, Xi 



and X2), and governed by the reactions [1 



A 
Xi 
Xi+B 
2X1 +X2 + C 



Xi + A, 
0, 

X2 + B, 

3X1 + c. 



(1) 

(2) 
(3) 
(4) 



These reactions conserve the numbers of molecules of 
types A, B and C in the system, while those of Xi and 
X2 are the dynamical degrees of freedom. The role of 
the substances A, B and C is mainly to set the rates 
with which the first, third and fourth reaction occur, re- 
spectively. 

The concentrations of the A and C molecules will be 
held constant in time in all variations of the model that 
we will consider, while the concentration of substance B 
will be used to apply an external driving force. The pre- 
cise manner in which this forcing is implemented will be 
detailed below. On the deterministic level, the system is 
described by the following two coupled ordinary differen- 
tial equations [H, [U 



±1 = 1 — Xl (1 -I- b{t) — CX1X2) ; 

±2 = X\ {b{t) - CX1X2) , 



(5) 



where xi{t) and X2{t) describe the time-dependent con- 
centrations of substances Xi and X2 respectively, the 
constant c the concentration of the C molecules (the 
concentration of the A molecules has been set equal to 
unity), and where b{t) is the externally controlled concen- 
tration of B-molecules. The unforced Brusselator is re- 
covered by setting b(t) = bo independent of time. In this 
unforced case the system may exhibit both fixed points 
and limit cycles, depending on the choice of the coeffi- 
cients bo and c (see 18] and references therein for details). 
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but no oscillatory approach to the limit cycles is possi- 
ble as discussed below. For later convenience we rewrite 
Eqs. ([5]) as X = A(x, t) where 



Ai{x,t) = 1 ~ xi {1 + b{t) ~ CX1X2) , 
A2{x,t) = X I {b{t) - cx 1X2) ■ 



(6) 



To complete the definition of the model it remains 
to specify the functional form of the forcing. We will 
here use a deterministic, periodically varying forcing, 
b{t) = 60(1 +ecos{ilt)), in Eqs. ([5]). The non-negative 
control parameter e sets the amplitude of the external 
driving, and f2 is its angular frequency. We restrict our- 
selves to e < 1 so that the concentration of i?-molecules 
remains non-negative. For e = we recover the unforced 
Brusselator. 



B. Deterministic dynamics: Floquet analysis and 
stability of limit cycles 

1. Floquet theory 

For periodic functions b{t + Tq) = b{t), and assuming 
£ 7^ 0, any periodic solutions of Eqs. ([5]) must have a 
time period T = nTa, where Tn = 27r/r2 and n is a 
positive integer. Numerical integration indeed shows that 
such cycles are found, though not for all choices of the 
model parameters. Furthermore, for the parameters that 
we have tested we only find limit cycles corresponding 
to n = 1. The stability of these periodic solutions may 
then be analyzed within the framework of Floquet theory 
P, . Assuming model parameters are set such that a 
periodic solution, x(t), of Eqs. ^ exists, one considers 
a small perturbation, ^{t), about this solution. To linear 
order one then has 



dm 
dt 



mm. 



(7) 



where K{t) denotes a 2 x 2 matrix with entries 
Kij{x.{t),t) = djAi{x.{t),t), i,j = 1,2 and dj denotes a 
derivative with respect to Xj . The explicit form of K(t) is 
given by Eq. (|Al|l of Appendix El but given that b{t) and 
x{t) are of period Tn it follows that K{t + Tq) = Kit). 
Equation (|A1[) is identical to that for the unforced case 
[l8| . except that here b is replaced by a time-dependent 
function b{t). 

In essence, Floquet theory states that, provided X{t) 
is a fundamental matrix of ([7]), then there exists a non- 
singular constant matrix B such that 



X{t + T,^=X{t)B, 



(8) 



for all t. While the Floquet matrix B will, in general, 
depend on the choice of the particular fundamental ma- 
trix X{t)^ its eigenvalues can be shown to be independent 



of this choice [24I . The eigenvalues of B are usually re- 



For the forced Brusselator model there are two multipli- 
ers, and we denote them by pi and p2 in the following. 
Since B is real, if one of the multipliers is real, so is the 
other. This is the situation found in two-dimensional au- 
tonomous systems. Characteristic exponents pi and /i2 
are then defined by pi = e'^'-^" for i = 1,2. 

We will now briefly discuss the properties of the result- 
ing Floquet multipliers. In order to make contact with 
the unforced system it is useful to distinguish between 
the cases 60 < 1 + c and 60 > 1 + c, as the attractor of 
the unforced system is a stable fixed point in the former 
case, and a limit cycle in the latter [l8|. 



2. The case 60 < 1 -I- c 

A trivial application of Floquet theory is to the un- 
forced case (e ^ 0), so that b{t) = 60. For 60 < 1 + c 
the deterministic system is then known to approach a 
fixed point, see e.g. [3 for further details. Floquet the- 
ory remains formally applicable as the matrix K{t) in 
Eq. ([7]) becomes time-independent at the fixed point; we 
will write K{t) = K* . Indeed, in this case, formally the 
time period of the matrix K can be set arbitrarily, as 
one has K{t + t) = K{t) for all t and t. Solutions to 
([7]) may be obtained directly by integration, and they 
can be written as ^(t) = exp{isr*t}^Oj where we have set 
the initial condition to be ^(0) = ^o- Considering two 
solutions, generated from two linearly independent ini- 
tial conditions, we can construct a fundamental matrix, 
X{t). It then follows from the form of the solutions to 
Eq. ([7]), and from Eq. ([5]), that the Floquet matrix B 
depends on the choice of the period, r, as 



B{r) 



(9) 



Denoting the eigenvalues of K* by A;, i = 1, 2, and those 
of B{t) by jOi(T), Eq. © yields the relation pi{T) = 
exp{AiT}. If the eigenvalues \i are complex, then they 
are a complex conjugate pair, A±. Setting c = 1 (which 
we do from now on), the eigenvalues of K* are given by 
A± = (foo/2) — 1 ± i\/6o(4 — 60), i-e. they are a pair of 
complex conjugates with non-zero imaginary part so long 
as 60 < 4. The imaginary parts of \± will be denoted by 
±w*. We will refer to uj* as the natural frequency of the 
unforced Brusselator. When forcing is applied, then in 
the limit e ^ 0, the functions Pi{T) = p±{t) are logarith- 
mic spirals in the complex plane, for bg < 1 + c, i.e. they 
are of the form 



p±(r) = e^(b„/2)-i]r (cos(tj*T) ± isin(w*r)) , 



(10) 



ferred to as the Floquet multipliers of the system ([7]). 



Following Wiesenfeld [l3|, we illustrate the position of 
the Floquet multipliers in the complex plane on an Ar- 
gand diagram, see Fig. [TJ The dashed line here cor- 
responds to Eq. (fTO|) at bo — 1.8 for the range r S 
[(7r/^*),(27r/^*)]. 

If the forcing amplitude e is small, but non-zero, the 
deterministic dynamics ([5]) no longer approaches a fixed 
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FIG. 1: (Color online) An Argand diagram of the Floquet 
multipliers, pi, for 6o = 1-8 and c = 1. The dashed red line 
indicates the position of the multipliers as e — > 0, with r € 
[tt/o;*, 27r/tJ*] (cf. Eq. (|10p ). The blue dot-dashed line shows 
the location of the multipliers at SI = 1.3, where different 
points on the line correspond to different values of e > 0. 
The solid line is for Q. = 1.7. The forcing amplitude is varied 
in the range e £ [10-^ 1]. The shaded area is the unit disk. 



point, but instead it is found to have a limit cycle. In this 
limit however, the deterministic trajectory is observed to 
remain close to the fixed point of the unforced case. The 
matrix K(t) in Eq. ([7]) then approaches K* as e — > 0. It 
then follows that pi expjAiTs^}, so that the Floquet 
exponents fii ^ Xi as e 0. We find from a numerical 
integration of the deterministic dynamics that increas- 
ing the level of the forcing amplitude tends to make the 
forced limit cycle more stable; that is, we find that the 
modulus of pi and p2 decreases when e is increased, as 
shown in Fig. [TJ 

Let us end this subsection by returning to the inter- 
pretation of complex Floquet multipliers. According to 
Floquet theory, a solution to Eq. ([7]) may be written as 
a linear combination of solutions which have the prop- 
erty ^i{t + Tu) = Pi^i{t) for i = 1,2. When the pi are 
complex conjugate pairs, this means that linear displace- 
ments from the periodic solution x(i) return to the limit 
cycle in elliptical spirals, in a way similar to the stable 
fixed point of the unforced case. We illustrate this typical 
behavior of complex Floquet multipliers in Fig. [21 



3. The case &o > 1 + c 

The case in which 6o > 1 + c is slightly more compli- 
cated than the one in which the unforced deterministic 
system approaches a fixed point. For 6o > 1 + c the un- 
forced system has a stable limit cycle solution 13 ; we will 
denote its angular frequency by loq, where ujq generally 
depends on bo and on c. One of the Floquet multipliers 




FIG. 2: (Color online) The top panel is a schematic plot of a 
deterministic approach, shown as a thin gray (red) curve, to 
a limit cycle with complex Floquet multipliers; the cycle itself 
appearing as the closed dark curve (blue). The lower panel, 
showing a stroboscopic section, illustrates the spiraling return 
to the limit cycle with a frequency distinct from that of the 
limit cycle itself. 

is equal to unity [3, [2l|, pi = 1, while the other one 
is found to be in the range < p2 < 1, consistent with 
a stable limit cycle attractor. We were not able to find 
any stable periodic solutions when integrating Eqs. ([5]) 
at small, but non-zero, forcing amplitudes e at generic 
forcing frequencies. At fixed values of bo and c, periodic 
solutions are however found for all fl when the forcing 
amplitude exceeds a critical value, which we denote by 
EciS^), suppressing a potential dependence on bo and c. 
For e > ed^) these solutions are stable limit cycles, and 
the corresponding Floquet multipliers lie within the unit 
circle. Here we will exclusively focus on this regime. At 
e — ed^) the multipliers have a modulus of one, so that 
the cycle loses its stability, and as in the previous subsec- 
tion, increasing the forcing amplitude reduces the moduli 
of pi and p2, as shown in Fig. [3] For our purposes it is 
sufficient to go on to study the case where the Floquet 
multipliers remain inside the unit circle, and to analyze 
the power spectra of stochastic fluctuations about the 
limit cycle in this regime. 



C. Stochastic dynamics and system-size expansion 

1. Specification of the Model 

We now turn to a discussion of the stochastic micro- 
scopic Brusselator system, as defined by the reactions 
([II-Q. Labeling the reactions by ly — 1, ... ,4, we de- 
note the rates with which each of the reactions occur by 
Ti^(n,t). These rates depend on the state of the system 
n = (77,1,712), where rii is the number of molecules of 
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FIG. 3: (Color online) Location of the Floquet multipliers in 
the complex plane for bo = 2.1 and c = 1. The dot-dashed 
blue line is for a forcing frequency of fl = 1.3, the solid black 
line for Q = 1.7. Periodic solutions are found above £c = 0.068 
and Ec = 0.14 respectively. Floquet multipliers are shown for 
Ec < e < 1 for both values of Q. The shaded area is the 
unit disk, Floquet multipliers approach the unit circle as e 
approaches Sc from above. 



species Xi, and for the forced system have an additional 
explicit dependence on time. For the Brusselator system 
Ti = N, T2(n) = ni, T3(n,t) = bo{l + ecos(m))ni and 
T4(n) = cN^^n\n2- The combinatorial factors are as in 
the unforced case [ll|. The time-dependent expression 
for T3 reflects the periodic forcing, implemented through 
an externally-controlled variation of the number of B- 
molecules in the system. We also define the vectors v^, 
u = 1, . . . , 4, each capturing the effects of a single occur- 
rence of a reaction of type v on the numbers of Xi and X2 
molecules in the system. For the Brusselator Vi = (1,0), 
V2 = (-1,0), V3 = (-1, 1) and V4 = (1,-1) 



2. Analytical description and system-size expansion 

The time evolution of the probability, Pn{t), of finding 
the system in state n at time t is then governed by the 
master equation 



4 

E 



[T,(n - v„ i)P„-v„ (t) - nin, t)Pn{t)] . 

(11) 

Solving the master equation analytically is generally 
not feasible, but an effective description in terms of a 
Langevin equation, valid at large, but finite, system size 
can be obtained by means of a van Kampen expansion in 
the inverse system size [23i] . 



This procedure is well-established and has been applied 
to a number of microscopic interacting particle systems, 
so that we do not describe the mathematical details here, 
but instead refer to [H, [2^ . The main idea is to expand 
realizations n(i) of the microscopic dynamics about a 
deterministic trajectory, x(t). 



^.5.(t) + -^|(t), 



(12) 



and to derive an equation of motion for the fiuctua- 
tions, ^{t), from an expansion of the master equation 
(fTTj) in powers of N~^/'^. To lowest order one finds that 
self-consistency requires x = A(x, i), where A(x, i) = 
(Ai(x, i), A2(x, t)) is given by the expressions in Eq. ([6]), 
recovering the deterministic dynamics of Eqs. ([5]). These 
equations may also be derived by defining 



(n(i))=5]nP„(0, 



and noting that 

d(n(t)) 
dt 



^v,r4(n(t)),i). 



(13) 



(14) 



where we have used a deterministic approximation to 
write (Tj.(n, i)) = T^{{n{t)) ,t). Equations ^ are then 
recovered by setting x(t) = {n{t)) /N. At next-to-leading 
order the van Kampen expansion gives a linear Langevin 
equation for the fluctuations, ^{t), about the determinis- 
tic trajectory which has the general form [l^, [2^ 



dm 
dt 



= K{t)m + Ht), 



(15) 



where, for the forced Brusselator, the matrix K{t) is de- 
fined in Eq. (|Aip . The term i{t) on the right-hand side 
represents a Gaussian noise of zero mean and with cor- 
relator 



(16) 



The matrix D{t) may be straightforwardly calculated 
from the van Kampen expansion p^ . [2^ . The explicit 
form for the forced Brusselator is given by Eqs. (|A2[) and 
(|X3)) in Appendix El 

Equation (jl5[) is a linear Langevin equation, and an- 
alytical progress is therefore possible. Of particular in- 
terest to us here are the correlation functions and power 
spectra of the fluctuations ^{t). The time-averaged ele- 
ments of the covariance matrix Cij{t,t') = {£,i{t)£,j{t')) 
are defined as 

C^Ar)^7^ r dt Um{t + T)). (17) 

We will in the following mostly focus on the diagonal 
elements Cii{T). Even though Eq. (|15p is linear, the an- 
alytical computation of Cu^t) requires several interme- 
diate steps, and final expressions need to be evaluated 
numerically. The details are left until the general theory, 
applicable to systems in an arbitrary number of dimen- 
sions, is explained in Sec. IIVI 



3. Comparison with simulations 

In Fig. [4] we compare results from the analytical calcu- 
lation just described, with measurements obtained from 
simulations of the microscopic dynamics. Simulations are 
carried out using the Gillespie algorithm suitably 
modified to account for the explicit time-dependence of 
the reaction rates induced by the external forcing [2^ . 
Measurements in simulations are taken after a suitable 
equilibration period in order to minimize the effects of 
transients. Fig. |4] shows results from the theory (lines) 
and from simulations (markers) , and as seen in the figure, 
the agreement between them is excellent, both for the 
time-averaged autocorrelation functions Cii(r) and the 
corresponding power spectra. The latter are obtained as 
the Fourier transforms of the correlation functions: 



P,(w) = J dre'^'CuiT). 



(18) 



In the numerical simulations we first measure Cn{t, t'), 
and then perform a time- average to obtain Cii{T). Sub- 
sequently a discrete Fourier transform is taken, to ob- 
tain Pi{ijj). From a practical point of view, Cii{T) is 
found only for r > 0, and then the even nature of the 
function (discussed later) invoked. Wiesenfeld sug- 
gested peaks would be expected to be seen at frequen- 
cies nfl ± Im/i, where n is a positive integer and Im/i 
is |Im/ii.2|, where /ii.2 are the two Floquet exponents. 
However, our results indicate that the presence or oth- 
erwise of such peaks depends strongly on the choice of 
model parameters, and in particular on the position of 
the Floquet multipliers in the complex plane. For the 
case shown in Fig. [H for example, pi,2 = 0.023 ± 0.46« 
and marked peaks are found at nfl — Im/x, but not at 
nO -|- Im/z. A second example is shown in Fig. [51 where 
we show data for a number of model parameters, result- 
ing in Floquet multipliers much closer to the unit circle 
than for the example shown in Fig. [?1 Peaks are now 
found at all nQ ± Imfi, with the peaks becoming more 
pronounced as the Floquet multipliers approach the unit 
circle (from within). In the limit |pi,2| 1, the relax- 
ation of autocorrelation functions becomes very slow and 
so larger values of r need to be taken into account when 
performing the Fourier transform. This makes both the 
analytical expressions and the Gillespie simulations more 
computationally expensive and, for the parameters illus- 
trated in Fig. [5l Gillespie simulation is not feasible. 





FIG. 4: (Color online) The autocorrelation (top panel) and 
the power spectra (bottom panel) of stochastic fluctuations 
about the deterministic trajectory of the forced Brusselator. 
Simulation results for fluctuations of the number of Xi- 
molecules are shown as open circles, while those for ^2 are 
indicated by full squares. The solid lines are the predictions 
of the theory, and are seen to match the simulations per- 
fectly. Model parameters are bo = 2.1, c = 1,0 = 1.3 and 
£ — 0.14. The corresponding non-trivial Floquet multipliers 
are p = 0.023 ± 0.46j, so that \p\ = 0.46. The imaginary part 
of the Floquet exponent is Im p = 0.32. The system size in 
simulations is A'' = 2 x 10^ and averages over 5000 indepen- 
dent runs are taken. Vertical lines in the lower panel mark the 
frequencies of nO, — Imp (solid lines) and nO, + Imp (dashed) , 
where n is a positive integer. 



III. WILLAMOWSKI-ROSSLER SYSTEM 



A. Microscopic model 

We have seen that forcing the two-dimensional Brus- 
selator opens up the possibility of complex Floquet mul- 
tipliers. This was not possible in the unforced case since 
there the deterministic dynamics is autonomous, leading 
directly to a Floquet multiplier of unity. Therefore, in 



order to see the effects of complex Floquet multipliers in 
an autonomous system, the simplest case has three di- 
mensions. One such system is the Willamowski-Rossler 
model and we shall study the particular form given in 
[2^ . [27| . The model may be written as a chemical re- 
action system, involving three species Xi, X2 and ^3, 
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Pico) 




FIG. 5: (Color online) Power spectrum of stochastic fluctua- 
tions, ^1, in the forced Brusselator system as obtained from 
the analytical calculations. Model parameters are bo = 2.1, 
c = 1, Q = 1.3, resulting in Sc = 0.068. The different curves 
correspond to forcing amplitude e = 0.07, 0.11, 0.15, from top 
to bottom, at the peaks. The corresponding Floquet multi- 
pliers have modulus 0.97, 0.66 and 0.41 respectively. Vertical 
lines are given at frequencies of nQ — Im^ (solid lines) and 
nCl + Im/i (dotted), where n is a positive integer and where 
Im/i ~ 0.31 for all three cases. 



defined by 







2X1, 


(19) 






2X2, 


(20) 






0, 


(21) 




1 


2X3, 


(22) 


^3 




0. 


(23) 



The parameters above and below the arrows indicate the 
relative rates with which each of the reactions occur. 
Absorbing potential combinatorial factors in the defini- 
tion of the model parameters fei , 62 , di , and d2 , the four 
reactions given in Eqs. (fT9|) and ((20|) occur with rates 
ri(n) = hm, T2{n) = dxn\N-\ Tsin) = &2n2 and 
74(11) = d2n\N^^. In isolation, these four reactions, (flQl) 
and (1201) , ensure that the average numbers of species Xi 
and X2, are of the order of iV, so that N is again a mea- 
sure of the system size. We will take the annihilation pro- 
cess (pTj) to occur with rate T5(n) = nin2N~^; the prcf- 
actor in the rate of the reaction is taken to equal unity in 
order to agree with [2^[23| . The mathematically interest- 
ing limit is that in which the number of X3 particles, 71-3, 
is of order N as well. This is the case when the remaining 



reaction rates are scaled suitably with iV. Specifically we 
will assume that occurs at rate T6(n) — nin^N^^ 
and at rate 17(11) — d^n^. The vectors, v^,, that 
correspond to the reactions i> = 1, . . . ,7 are given by 



VI -(1,0,0), 
V3 = (0,1,0), 
V5 - (-1,-1,0,) 
V7 = (0,0,-1). 



V2 = (-1,0,0), 
V4- (0,-1,0), 
V6 = (-1,0,1), 



(24) 



B. Deterministic Dynamics and Frenet Frame 

As in the case of the forced Brusselator, we may now 
find the equations of the corresponding deterministic dy- 
namics using Eq. For the Willamowski-Rossler 
model these are 

a;i=^i(x) = xi{bi - dixi - X2 - X3), (25) 

i2=A3(x) = X2{b2 - d2X2 - Xl), (26) 

i3=^3(x) = xsixi-ds). (27) 

There are a total of six fixed points of this system, but 
only one at which all concentrations are non-zero. This 
fixed point is given by 



X* = d3. 



- dids 



bo - rf3 



(28) 



The stability matrix, Kijlx) — djAi{'x.), at this fixed 
point may be found from Eqs. (|A4[) and (jASP in Ap- 
pendix 1^ by setting x.(t) — x* . If the above non-trivial 
fixed point is unstable then limit cycle solutions of the 
deterministic equations may exist. Such solutions have, 
for example, been reported in (26l. [27j. and we will focus 
on this limit-cycle regime in this section. 

Following the notation of the previous sections we will 
denote the deterministic limit cycle trajectory by x(i), 
and we will write ^(<) for the fluctuations about it, again 
as before. Much of the formalism we require has either 
been discussed in SecUll or in [ll|. In particular one has 
an equation of the form ([7]) within a linear stability analy- 
sis of the limit cycle, and in the absense of noise. A direct 
consequence of the system being autonomous is that the 
velocity vector, x(i), is itself a solution to Eq. ([7]). Since 
the velocity is periodic, one of the Floquet multipliers is 
equal to unity, as it is generally the case for limit cycles of 
autonomous systems. The dynamics is marginally stable 
in the direction of the velocity, so that longitudinal fluc- 
tuations behave diffusively and may grow without bound 
in the long run [3 HIl • We will focus our interest instead 
on the fluctuations in the transverse directions, since it is 
these that have the oscillatory behavior of interest to us. 
For stable limit cycles and in the absence of persistent 
noise, these transverse fluctuations decay in a manner 
characterized by the remaining Floquet multipliers. If 
the latter are complex, and if the system is subject to 
intrinsic noise, as induced by the underlying microscopic 



8 



dynamics at finite system sizes, we expect these fluctu- 
ations to be enhanced into quasi-cycles about the hmit 
cycle. 

In order to separate longitudinal from transverse 
modes we need to introduce a suitable frame of refer- 
ence. Such co-ordinates are provided by the Frenet frame 
[2^ . which may be constructed by applying the Gram- 
Schmidt orthogonalization procedure to the first three 
time derivatives of the limit cycle solution x(i). Specifi- 
cally, the co-moving basis vectors ei{t), i — 1, 2, 3 of the 
Frenet frame are constructed sequentially, as discussed 
further in Appendix |B] The fiuctuations are governed by 
a Langevin equation of the form (jlSp . In order to isolate 
the transverse fluctuations, we rotate the Langevin equa- 
tion into the Frenet frame. After the rotation, defined by 
a matrix, J(t), the Langevin equation takes the form 

ci{t)^K'°\t)qit)+git), (29) 

where we follow our earlier paper [Tsj and write q(t) = 
J{t)^{t) for the fiuctuations in the Frenet frame. The ma- 
trix is periodic and given by K*'°^{t) = J{t)K{t)J^^{t) + 
j{t)J-^{t) (see Appendix E]) and g{t) = J{t){{t) is the 
rotated noise term. It follows from Eq. (fT6)l that the com- 
ponents of g{t) are each Gaussian white noise variables 
with zero mean and correlators 

{gS)g,{t')) ^2G,,{t)S{t~t'), (30) 

where G{t) = J{t)D{t)J-^{t). 

For autonomous systems, it is shown in Appendix |B] 
that the existence of a longitudinal direction as described 
above, implies that the elements of the first column of 
the matrix K^°'^{t) vanish, except for the entry in the 
first row. A consequence of this is that the transverse 
dynamics may be effectively considered independently 
of the dynamics in the longitudinal direction. For the 
Willamowski-Rossler limit cycle this yields a pair of cou- 
pled linear Langevin equations in the two transverse di- 
rections, with exactly the same mathematical form as 
those of the forced Brusselator model. Hence the same 
techniques as before may be applied to produce analyti- 
cal curves for the autocorrelations and power spectra in 
the two transverse directions. Note that in our previous 
work (is] , we were able to simplify the rotated Langevin 
equations further by a rescaling of the coordinates in the 
Frenet frame. We do not apply this additional transfor- 
mation here, since for our purposes it is not essential. 

C. Stochastic Simulation and Results 

The Gillespie algorithm can again be used to gen- 
erate realizations of the microscopic dynamics defined 
by Eqs. (|19p - ([23|) . Since one Floquet multiplier in the 
Willamowski-Rossler system is equal to unity, there is a 
diffusive mode in the longitudinal direction. This means 
that the time-evolution {ni{t)/N,n2{t)/N,n^{t)/N) of 
any single realization of this stochastic process may not 



remain close to the deterministic trajectory x(i), but in- 
stead (|n(t)/A^ — x(i)p) ~ t, where | • | stands for the 
Euclidean norm. This complication is not present in the 
driven Brusselator discussed in Sec. [Tll since in that case 
no such longitudinal diffusive mode exists. 

This issue can however be dealt with as discussed in 
[l8l |. The procedure of extracting the deviation from 
the limit cycle is as follows: for every given data point 
n{t)/N generated by the Gillespie algorithm one iden- 
tifies the point x(n(t)) on the limit cycle trajectory 
which is geometrically closest to n{t)/N, and then uses 
K{t) — n{t)/N — x(n(<)) as the displacement vector. As 
described in [T^ the longitudinal component of K,{t) van- 
ishes, i.e. one has x.k = 0, while the remaining compo- 
nents define a stochastic process in the co-moving trans- 
verse plane, and as seen in [Tsl the magnitude of k re- 
mains of order 7V^^/^. This procedure allows one to ef- 
fectively decouple the diffusive longitudinal mode from 
the transverse ones, and we will focus on the transverse 
components in the following, in order to characterize 
stochastic oscillations about the deterministic limit cycle. 
These components are then expressed in the Frenet co- 
ordinates, defined at x(n(i)). As an illustration, trajec- 
tories of the transverse components obtained from a sin- 
gle realization of the microscopic dynamics are shown in 
Fig. [S] for a fixed set of model parameters. In this figure, 
N{t) denotes the normal component, N{t) — K{t).e2{t), 
and B{t) denotes the deviation from the limit cycle in 
the binormal direction, B{t) = K{t).e3{t). Recall here 
that §2 and 63 define a co-moving frame, i.e. that they 
carry a time-dependence as well. 

In Fig. [3 we show the resulting power spectra, and 
find very good agreement between simulation and the- 
ory for both the normal and binormal directions. There 
is a slight systematic deviation of data points from the 
theory, which occurs at integer multiples of the limit cy- 
cle frequency. We attribute these to remnants of the 
deterministic dynamics. The data shown in Fig. [7| was 
taken at model parameters resulting in complex Floquet 
multipliers with a modulus of approximately 0.3, and 
peaks are found in the power spectra close to frequencies 
nwoilni/Lt, where wq is the angular frequency of the limit 
cycle. However, we also note that peaks are not observed 
at all frequencies nujQ ± Im/i, especially in the spectrum 
of normal fluctuations. As for our findings in the driven 
Brusselator, this may be due to the fact that the Floquet 
multipliers in the example shown in Fig. [7] are relatively 
distant from the unit circle in the complex plane. Again 
based on our observations in the driven Brusselator one 
may expect additional peaks at frequencies nujQ ± Im^ to 
emerge as the Floquet multipliers move closer to the unit 
circle. Despite an extensive search we have however not 
been able to find a set of model parameters which would 
result in Floquet multipliers with modulus close to unity, 
so that we are not able to give any further confirmation 
of this expectation here. We conclude this section by re- 
iterating our main result, the near perfect agreement of 
the analytically obtained power spectra with simulations 
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FIG. 6: The fluctuations in the directions transverse to 
the limit cycle trajectory in the Willamowski-Rossler model. 
Data for the normal component, N(t), and the binormal di- 
rection, B{t), are shown for a single realization of the stochas- 
tic simulation. Model parameters are bi — 80, 62 = 20, 
di = 0.16, d2 = 0.13, and ds = 16. 



1000 - 



Pb(co) 




FIG. 7: (Color online) Comparison of the theoretical and 
simulated estimates for the power spectra of transverse fluc- 
tuations in the Willamowski-Rossler system. The top panel 
shows the normal fluctuations (in the direction 62) while the 
bottom panel compares those in the binormal direction, ea. 
Model parameters are again set to 61 = 80, &2 ~ 20, di — 0.16, 
d2 = 0.13, and da = 16. Vertical lines are given at frequencies 
nujo + Im jj, (dotted) and nujo — Im fi (solid), with n a positive 
integer. The numerical value of ujq is 17.25, and the non- 
trivial Floquet multipliers are p — —0.002 ± 0.303 (resulting 
in IpI = 0.30 and Im = 4.33). 



as shown in Fig. [T] 



IV. GENERALIZATION TO HIGHER 
DIMENSIONS AND THE COUPLED 
BRUSSELATOR 

A. General theory 

It is expected that in the study of any real system, 
for example in biochemistry or in ecology, the number 
of distinct species, S, would be significantly larger than 
two or three. It is also possible that solutions to the S*- 
dimensional deterministic equations in such a model may 
be periodic orbits, x(t). Hence, in this section we demon- 
strate the natural extension of the analysis in the previ- 
ous sections to models of arbitrary dimension. Whether 
or not the system is autonomous, we begin with the van 
Kampen system-size expansion which yields a set of S 
coupled and linear Langevin equations for stochastic fluc- 
tuations, ^{t). We simply note that their form naturally 
extends to arbitrary dimension and is unchanged from 
(fT5|) . where the S x S matrix K{t) for the drift is given 
by Kij{t) = Kij{x{t)) — dAi{x.{t))/dxj, and the sym- 
metric S X S matrix for diffusion, D{t) = D{x{t)), is 
calculated from the system-size expansion. 

The subsequent steps of the analysis then depend on 
whether the system under consideration is autonomous 
or not. For non- autonomous systems no rotation is re- 
quired, and one proceeds directly with the Langevin 
equation in Cartesian co-ordinates in S dimensions. 
If the system is autonomous, as in the case of the 
Willamowski-Rossler model, one first needs to rotate into 
the S'-dimensional Frenet frame, and then to separate 
off the longitudinal component, resulting in a Langevin 
equation in 5 — 1 dimensions for the transverse com- 
ponents. The Frenet frame is defined in S dimensions 
in Appendix |B1 This then specifies the rotation matrix 
= (ei, . . . ,65), which is evaluated on the limit cycle 
so that J{t) = J(x(t)). The formalism is a straightfor- 
ward generalization of that described in Sec. IIIll for the 
Willamowski-Rossler model, except that there are now 
{S — 1) transverse directions, rather than just two. 

Thus, for both autonomous and non-autonomous sys- 
tems one eventually ends up with a Langevin equation 
in d dimensions, where d = S — 1 for the autonomous 
case, and d = S for non-autonomous systems, such as 
the driven Brusselator. The further steps of the calcu- 
lation can hence be discussed simultaneously for the au- 
tonomous and non- autonomous cases. As described in 
more detail in Appendix [C| the solution of this Langevin 
equation can be expressed in terms of any fundamental 
matrix X{t) of the corresponding homogeneous equation. 
Since the drift matrix, K{t) (denoted by K(t) in the au- 
tonomous case) is periodic, Floquet theory |22| asserts 
that a canonical fundamental matrix may be written in 
the form X{t) = P{t)Y{t), where P{t) and Y{t) are d x d 
matrices. The matrix P{t) is periodic with the same pe- 
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riod as the drift matrix while the matrix Y{t) is given by 
Y{t) = e'^'^s{A'i,. .,Md}t^ where the fi,, i ^ I, . . . , d are the 
Floquet exponents of the d x d homogeneous system. 

The periodic matrix, P{t), acting on the left is, in ef- 
fect, a transformation matrix from the Floquet solutions 
to the coordinates of the Langevin equation, while its in- 
verse makes the reverse transformation. The matrix Y{t) 
is a diagonal exponential matrix with entries e'^'*, with 
Re/ii < 0, for all i, for a stable limit cycle. It therefore 
acts on different Floquet solutions in different ways, re- 
ducing the value of some more quickly than others. The 
general solution of the Langevin equation (jCip which we 
wish to analyze, can be given in terms of the matrices 
P{t) and Y{t) and is given explicitly by Eq. I|C5|) in Ap- 
pendix [Cl 

Given the symmetric and periodic noise matrix, D{t) 
in the non-autonomous case — which we generally denote 
by G{t) using the notation of the autonomous case — 
we may calculate the autocorrelation function in closed 
form. In the basis corresponding to the Floquet solu- 
tions G{t) becomes the symmetric and periodic matrix 

r{t) = P-^{t)G{t){p-^)'^ . These noise contributions 
are then integrated over one time period of the deter- 
ministic limit cycle, T, but weighted by decaying expo- 
nentials from the Y{t) matrix, to yield another symmet- 
ric and periodic matrix, A(t) (see Eq. (|C7[) ). which gives 
the various covariances of the fluctuations in the space 
of the Floquet solutions. However, the focus of our in- 
terest is in the two-time correlations of the fluctuations 
which are shown in Appendix [C] to equal C{t + T,t) = 
2P{t + T)Y{T)K{t)P^{t). Therefore the autocorrelation 
function itself equals 



the reactions, 



C{t) 



P{t + T)Y{T)A{t)P'^{t)dt. 



(31) 



for r > 0. The diagonal elements of C(t) turn out to 
be even functions of r, as they ought to be. Power spec- 
tra, Pi{uj) for i — l,...,d, may then be calculated as 
the Fourier transform of diagonal elements of C(t), as in 
Eq. nil). 



B. The case of two coupled Brusselator systems 



A 

ATa + X2 
2X3 + A4 + C" 



A3 + A, 

A4 + A2 , 
3 A3 + C. 



(32) 
(33) 
(34) 
(35) 



Given that substance A is part of both units, the sec- 
ondary Brusselator therefore has the same system size as 
the primary one. The deterministic dynamics is given by 

Xl = 1 — a;i(l -I- 6 — ca;ia;2), (36) 
X2 — xi{b~cxiX2), (37) 
X3 — 1 — a;3(l -I- a;2 — c'x3a;4), (38) 

±4 = X3{X2 - c'x3X4). (39) 

When 6 > 1 -|- c there is a limit cycle in the primary 
Brusselator; we will again denote its angular frequency 
by LUo- These oscillations of the primary Brusselator act 
as a periodic forcing on the second, and for all parame- 
ters studied here, the second Brusselator shows cycles at 
the above frequency luq. The two Brusselators together 
form a four-dimensional autonomous system. Hence, we 
will study the fluctuations of the large system-size dis- 
crete system which act transverse to the limit cycle. In 
this example then, we discuss the normal, 62, binormal 
63, and trinormal 64, directions. Once the periodic drift 
K(t) and diffusion D{t) matrices, given by Eqs. (|A8[) - 
(|A11|) in Appendix 13 are rotated into the Frenet frame, 
we then calculate power spectra of transverse fluctuations 
via Eq. (pij) . The results of this are presented in Fig. [5] 
for the model parameters b = 3.3, c = 2, and c' = 1. We 
find very good agreement between theory and simulation 
performed using the Gillespie algorithm. For these pa- 
rameters, one of the non-trivial Floquet multipliers, p2, 
is real and positive, while the remaining two, p±, take 
on complex conjugate values. However, these multipliers 
are not associated with any particular transverse direc- 
tion, as can be seen from the power spectra: in all three 
directions, 62, e3, and 64, peaks are found at frequencies 
equal to a multiple of luq, but also at those associated 
with the imaginary parts of the complex Floquet expo- 
nents, ojQ + Im fi± . While the general formalism we have 
developed in this section has been illustrated on the con- 
crete example of the coupled Brusselator, it should be 
clear that it can be applied quite generally to investigate 
the fluctuations about a limit cycle in S'-dimensions. 



In order to demonstrate the method on a concrete ex- 
ample, we will study a model composed of two coupled 
Brusselator systems. Two Brusselator units can be cou- 
pled in a number of different ways and here we construct 
the coupling in such as way as to draw parallels with the 
forced Brusselator discussed earlier. Chemical species Ai 
and A2 form a primary Brusselator through reactions, 
([T|)-(l4|), with constant populations of A, B and C. We 
now also introduce species A3, A4 and C", which follow 



V. CONCLUSIONS 

The phenomenon of stochastic amplification due to de- 
mographic, or intrinsic, noise has been qualitatively un- 
derstood for fifty years, but it is only recently that it 
has been comprehensively and quantitatively described. 
This has been due in large part to the application of the 
technique of the system-size expansion, which is able to 
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FIG. 8: (Color online) Power spectra of the fluctuations about 
the limit cycles of the coupled Brusselator system. Data is 
shown for the three transverse directions, with black circles 
indicating the normal direction, dark gray (red) squares the 
binormal component, and light gray (green) diamonds the 
trinormal directions. Markers show results from simulations 
at a system size of 2 x 10^ and averaged over 2000 runs. The 
solid lines are from the theory, and as seen in the figure the 
agreement with simulations is near perfect. Vertical lines are 
shown at frequencies ncuo (solid), nujo +Im/i (dotted) and 
nujo — Im/i (dashed), with n a positive integer. 

reproduce results obtained by numerical simulations to a 
remarkable precision. In fact, although this method al- 
lows for a systematic expansion in powers of 1/ VlV, there 
is usually no need to go beyond next-to-leading order. In 
essence, application of the method means that the use of 
numerical simulations to understand the cycles induced 
by noise could be dispensed with entirely. 

If the systems under study are subject to an external 
periodic driving, for example biological systems subject 
to an annual cycle, then the deterministic dynamics may 
have a limit cycle as its stable state. In this paper we 
have investigated the effect that demographic stochastic- 
ity will have on this state. On general grounds one might 
expect that if the limit cycle was approached in an os- 
cillatory manner, then stochastic cycles about the limit 
cycle could be sustained. We have shown that once again 
the system-size expansion may be applied to gain a quan- 
titative understanding of this phenomena. The analysis 
is considerably more elaborate than in the case where 
the deterministic dynamics approaches a fixed point, but 
once again the method gives excellent agreement with 
numerical simulations. 

The signature for the oscillatory approach to limit cy- 
cles is that the associated Floquet multiplier should be 



complex. This can occur for nonautonomous systems in 
two or more dimensions or autonomous systems in three 
or more dimensions. Since the eigenvalues of a typical 
real matrix in these dimensions will generically be com- 
plex, one might expect complex Floquet exponents to 
be common. Our investigations of various models, al- 
though far from comprehensive, suggests that they are 
quite common in periodically driven systems, but not so 
common in autonomous systems that are generally stud- 
ied. There may be a dynamical reason for this, but it 
is as likely that this is due to the nonlinear systems ap- 
pearing in the literature being selected for their period 
doubling transition to chaos, rather for the structure of 
their limit cycles. 

In the past it was said that intrinsic noise could turn 
oscillatory decay to a fixed point into sustained oscilla- 
tions. It was expected that these oscillations would have 
periods Im , where were the eigenvalues of the sta- 
bility matrix for that fixed point. This is only true in a 
very broad sense, as studies over the last few years have 
shown. In reality the period may significantly deviate 
from Im Ai due to other factors, and the amplitude of the 
fluctuations may be much larger than might be expected 
due to a resonance effect. Analogously, one might guess 
that intrinsic noise could turn oscillatory decay to a limit 
cycle into sustained oscillations about that cycle and that 
these oscillations would have periods hloq ± Im fj.^ , where 
luq is the period of the limit cycle and /i^ are the Flo- 
quet exponents associated with that limit cycle. We have 
shown in this paper that this is indeed the case in a broad 
sense, but as for the case of the fixed point there is much 
more to the story than this. For instance, the expressions 
ncuQ ± Im fii are again just an approximation to the fre- 
quencies and the amplitude of the oscillations will vary 
significantly depending on a number of factors, such as 
the magnitude of the Floquet multipliers. Fortunately, 
the system-size expansion once again gives results which 
are in excellent agreement with simulations and gives us 
a way of exploring the nature of these fluctuations. We 
expect that the ideas presented in this paper will have 
a number of applications, which we hope to explore and 
report on in the future. 
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become periodic functions of time. For the remainder of 
this appendix, we shall suppress the time dependence of 
x(t) for greater clarity. 



Forced Brusselator 



[1 + b{t) - 2CX1X2] cxl 
\b{t) — 2CX1X2] 



D{t) = 



Di{t) -£>2(i) 



where 



= -{l + Xi[l + b{t) + CXiX2]}.. 
D2{t) = ^{xi[b{t) + CXiX2]}. 



2. Willamowski-Rossler Model 



where 



and 



-Xl -xi 



\ 



-X2 K22{t) 

\ X3 Kis{t)J 



Kii{t) = 61 - 2dixi -X2- X3, 

K22 (t) = 62 - 2d2X2 - Xl , 

Kziit) = xi-ds, 



m = 



/ Dn{t) Di2{t) Di3{t)\ 

Di2{t) D22{t) 
\Dis{t) D33it)J 



(Al) 



(A2) 



(A3) 



(A4) 



(A5) 



(A6) 



where 



APPENDIX A: EXPLICIT FORMS OF 
MATRICES 

The matrices which appear in the description of the 
fluctuations about the deterministic trajectory are given 
in this appendix. The drift matrix K{x.) and the diffusion 
matrix D(pi) are naturally functions of the concentration 
X. However, when the solutions of the deterministic dy- 
namics, x{t + T) = x{t), are limit cycles they themselves 



Dii{t) = -xi{bi+dixi+X2 + X3), 



Di2{t) 
Dl3{t) 



= -^XlX2, 



= --jXlXi, 



D22{t) = -X2{b2 + d2X2+Xi), 

Daait) = ^X3{xi + ds). 



(A7) 
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K{t) = 



where 



and 



3. Coupled Brusselators 

/ Ki{t) -1 cxf 
-Ki{t) -cxl 

-x^i K3{t)-1 , 

V X3 -K3{t) 

Ki (t) = 2cxiX2 — b, 
K^{t) = 2c'xzXi-X2, 

I D^it) -D2{t) 
-D2{t) D2{t) 







'C'xl 



(A8) 



(A9) 







V 











D:,{t) -Di{t) 
-Di{t) Di{t) ) 



(AlO) 



where in addition to (|A3[) we have, 
1 



[1 + 0:3 (1 + a;2 + c'x3a:4) 



^4(t) = -j^X3 {X2 + C X^Xi) . 



(All) 



APPENDIX B: THE FRENET FRAME 

In this appendix we wiU discuss the background, and 
develop the formahsm, relating to the co-moving frame 
which we use to study the fluctuations from the limit cy- 
cles discussed in the main text. Such a frame, called a 
Frenet frame [29| , is a natural way to study displacements 
from a deterministic trajectory in any number of dimen- 
sions. Here we will denote the number of dimensions by 
S. 

Consider the general autonomous problem which we 
may describe by a system of non-linear, homogeneous, 
first-order equations, 



dx 



= A(x). 



(Bl) 



Following |29|, we may define the Frenet frame by ap- 
plying the Gram-Schmidt orthogonalization procedure to 
the time-derivatives of the solution, x(t). So long as the 
time derivatives are linearly independent, this gives the 
basis vectors, e.i{t), of the frame to be 

d*x(<) ^ /d'xm „ , A . , , 



(B3) 



We may now construct the matrix which transforms 
from Cartesian co-ordinates to the Frenet frame to be 
J{t) = (ei(t), . . . , es{t))^ . The transformation is by con- 
struction an orthogonal matrix O(S'), and as such has the 
property that {t) = J~^{t) for all times. 

We now wish to consider the effect of this transforma- 
tion on the equation of a linear fluctuation, $,{t), about 
the deterministic solution, x(t). For the time being we 
will neglect the noise term and consider the homogeneous 
equation, ^(t) = K{t)^{t). The transformation to the 
Frenet frame takes the form ^{t) i~> q(t) = J{t)^{t). 
Then ^(t) = (i(t) -I- J(t)K{t))^{t) and so the rotated 
displacement obeys the linear equation. 



q{t) = K'°'{t)q{t), 
where K*^°\t) = K'{t) + R{t) and where 



(B4) 



K'{t) = J{t)K{t)J-\t), R{t) = j{t)J~\t). (B5) 

We now evaluate the elements of the first column of 
the matrix K^°^. These have an especially simple form, 
with Kl°^ = for j > 1. This follows from the fact 
that, for an autonomous system, x.{t) — K{t)x{t), and 
so the "velocity" x(t), is a solution of the homogeneous 
equation that we are considering. From this, and from 
ei(t) = ±{t)/\ic{t)\, it follows that 



1 



|x(i)l 



(B6) 



The second term in the definition of K*""^ (t) , R{t) , may 
be written in terms of the basis vectors and, due to their 
orthogonality properties, we have 

de,(t) dei(t) 



dt 



for i 7^ 1. The rate of change of the longitudinal basis 
vector is given by (x(t) — ei(ei • x(t))) /|x(i)| and so 

R^l{t) = -^efc(i) .#(0, i ^ 1. (B8) 
|x(i)| 



Adding Eqs. (|B6P and (jBSp . and noting that Rn = 0, we 
have 

Kl°\t) = > 1); Klfit) = -zi^m ■ m- (B9) 

So all of the elements of the first column of K*""^ (t) vanish, 
apart from the element which is also in the first row. The 
significance of this is that the transverse displacements, 
which we denote by r(i) decouple from the longitudinal 
displacements, denoted by s{t). So writing a general dis- 
placement as q(t) = {s{t),r{t)), we have 



r(t) = Kit)r{t), 



-Ksr{t)-r{t), (BIO) 
(Bll) 
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where the vector K.sr{t) is the (5— l)-dimensional vector 
Kl°^{t) and where K{t) now describes the purely trans- 
verse drift behavior. So the Frenet frame always sepa- 
rates the equation of motion for the linear fluctuations 
into longitudinal and transverse parts and the transverse 
motion is free from any influence by the longitudinal mo- 
tion. 



APPENDIX C: AUTOCORRELATIONS OF 
PERIODIC LANGEVIN EQUATIONS 

The equations which describe small perturbations 
about the limit cycle either have the form ([7]) for non- 
autonomous (forced) systems or the form (jBlip for au- 
tonomous (unforced) systems. In the latter case longi- 
tudinal displacements have been excluded, but once this 
has been done, the analysis for both cases is identical. 
So we can develop the theory for both together, we will 
adopt the notation of the autonomous case, that is, start 
from the equation r{t) = K{t)r{t). It should then be un- 
derstood that in the non-autonomous case the replace- 
ments r{t) ^{t) and K{t) K{t) should be made. 

The results of Floquet theory (2^ tell us that, when 
K{t + T) = K{t) for all t, one may generally flnd d lin- 
early independent solutions to the homogeneous equation 
Y{t) — K{t)r{t) which have the form r{t) — pi(i)e'^'*. 
Here /i^, i = 1, . . . ,d, are the Floquet exponents, which 
may in general be complex, and the functions Pi{t) are 
periodic with the period, T. From these solutions, the 
the canonical fundamental matrix, X(t), may be con- 
structed. It has the special property that the constant 
Floquet matrix, B = X^^{t)X{t + T), is diagonal with 
elements equal to the Floquet multipliers. Grimshaw 
[131 appends a subscript to denote the canonical choice 
which results in a diagonal Floquet matrix, but since 
we will only deal with such a choice in this paper, we 
omit this subscript. However when carrying out nu- 
merical work, it should be recognized that in general 
the solutions which are found will be linear combina- 
tions of solutions of the form pi(i)e'^'*. These can be 
used to flnd a (non-diagonal) B, the eigenvectors of 
which can be used to construct a similarity transfor- 
mation to a canonical form. An alternative way of de- 
scribing the canonical solutions is to deflne the periodic 
matrix P{t) = {pi{t), . . . ,Pd{t)) and the diagonal expo- 
nential matrix Y{t) — exp{diag(/ii . . . /i(j)0- terms 
of these the canonical fundamental matrix is given by 
X{t) = P{t)Y{t). 

Moving on to the fluctuations about the periodic so- 
lutions of the deterministic dynamics, the linear stochas- 
tic fluctuations obey a Langevin equation (jlSp . with the 
noise correlator given by Eq. (|16[) . for non- autonomous 
(forced) systems and a Langevin equation (|29p . with 
the noise correlator given by Eq. where G{t) = 

J{t)D{t)J^^{t), for autonomous (unforced) systems. To 
separate out the latter into longitudinal and transverse 



components, we note that in Appendix [B] we wrote 
q(i) — {s{t),r{t)), and now we analogously write g{t) — 
{gs{t),gr{t)). Then, since the transverse fluctuations de- 
couple from the longitudinal fluctuations, the Langevin 
equation for purely transverse fluctuations r{t) may be 
written as 



(CI) 



The noise correlator ([50)) can be expressed in terms of 
transverse and longitudinal components by decomposing 
G{t) as follows: 



Git) 



Gssit) Gj(i) 
G,,(t) G{t) 



(C2) 



Since the vector Gsr{t) is typically non-zero, the ran- 
dom variables, gs and g^, generally remain statistically 
correlated in the rotated frame. However, this is only im- 
portant if we intend to evaluate simultaneous values of 
both gs{t) and gr{t) and this we do not do, because we 
have already shown for the noiseless case that the trans- 
verse displacements are independent of longitudinal one. 
Therefore the only noise correlator we require is 



{gr{t)-gJ{t'))=2Gmt~t'). 



(C3) 



Once again we will develop the theory using the notation 
of Eqs. ((CT|) and ((C3)) . but it apphes equally to Eqs. 
and (fTCjl. 

Floquet theory may be applied to linear inhomoge- 
neous equations of the form (jCip . as well as to homo- 
geneous equations such as r{t) = K{t)r{t) [IJ. To solve 
Eq. (jCip . we proceed in the standard way and add a 
particular solution of the equation to a general solution 
of the corresponding homogeneous equation. This yields 



rit) 



X{t)ro+Xit) f X-\s)gr{s)ds, 



(C4) 



for t > to and with the initial condition r(to) — X{to)ro- 
Since we will not be interested in the effects of transients 
in this paper, we set the initial conditions in the inflnitely 
distant past, to ^ —oo. A change of integration variable 
s — > s' = t — s in the solution (|C4p now gives 



r{t) = P{t) / Y{s')p-\t - s')gr{t - s')ds', 



(C5) 



where we have used the fact that, since Y{t) is a diagonal 
exponential matrix, Y{ti + ^2) = ^(^1)^(^2)- 

Of course, r{t) is a stochastic variable, and we will 
typically be interested in flnding correlation functions, 
principally the two-time correlation function C(t + T, i) — 
{r{t + T)r^{t)). Taking r > 0, the solution jCKj) gives 



Cit + T,t) = 2P{t + T)Y{T)A{t)P^{t), 



(C6) 
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where we have introduced the symmetric and periodic 
matrix integral, 



A(i) = / Y{s)T{t- s)Y{s)ds, 
Jo 

and, in turn, the symmetric and periodic matrix 



T{s) ^ p-\s)Gis) {P~\s))- 



(C7) 



(C8) 



All of the functions in Eq. (|C6p are deterministic and 
may be evaluated given a good numerical estimate for 
the limit cycle solution x(t) . 

The infinite integral for A{t) may be evaluated as a 
re-summed finite integral due to the periodicity of r(s). 
The result, in terms of Floquet multipliers, pi, is then. 



1 



1 - PiPj Jo 



(C9) 



for i,j = 1, . . . jd. The origin of the prefactor is from 
an infinite geometric summation, X]n=o(^'*/^i)"' which is 
convergent when the Floquet multipliers are inside the 
unit circle. 



Finally, although the details are not presented here, 
an expression can be found for r < 0. It turns out that 
C(r) = C{—t)'^, as it ought. Hence the final form is 
given by Eq. (|3ip for r > 0, and can be found from 
Eq. ([3T|l for T < 0, supplemented by the condition C(t) — 
C(-r)T. 



